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The uhimate apphcation of Einstein's field equations is to empirically determine the geometry 
of the Universe from its matter content, rather than simply assuming the Universe can be repre- 
sented by a homogeneous model on all scales. Choosing an LTB model as the most convenient 
inhomogeneous model for the early stages of development, a data reduction procedure was recently 
validated using perfect test data. Here we simulate observational uncertainties and improve the 
previous numerical scheme to ensure that it will be usable with real data as soon as observational 
surveys are sufficiently deep and complete. Two regions require special treatment — the origin and 
the maximum in the areal radius. To minimize numerical errors near the origin, we use an LTB 
series expansion to provide the initial values for integrating the differential equations. We also use 
an improved method to match the numerical integration to the series expansion that bridges the 
region near the maximum in the areal radius. Because the mass enclosed within the maximum obeys 
a specific relationship, we show that it is possible to correct for a fixed systematic error in either 
the distance scale or the redshift-space mass density, such that the integrated values are consistent 
with the data at the maximum. 

PACS numbers: 04.20.-q, 95.30.Sf, 98.65.Dx, 98.80.Jk 



I. INTRODUCTION 

Einstein's equations relate the geometry of spacetime to the distribution of matter, and the most important ap- 
plication of this key physical concept is the study of the Universe. As the amount and precision of cosmological 
data improve, we will be able to determine not only the overall curvature of the Universe, but its detailed geometric 
structure too. 

It is generally assumed that the Universe is homogeneous above the scale of supcrclusters and that it can be 
represented by a Friedmann-Lemaitre-Robertson- Walker (FLRW) model. This assumption has served cosmology well 
and led to a very good understanding of many features of the observed Universe, which implies it cannot be seriously 
wrong. Yet it must be acknowledged that homogeneity is indeed an assumption that should eventually be checked 
against observations, once the data are sufRciently accurate and complete over a large enough range of redshifts. 
Because of our fixed position in the Universe, there are two distinct aspects to homogeneity — isotropy about our 
position, and uniformity with distance from us. Isotropy is relatively easy to check, but radial homogeneity has 
many intrinsic complications. Therefore a thorough verification will require considerable effort. The assumption of 
homogeneity is so pervasive, and underlies so much theoretical and observational work, that there is a real danger of 
a circular argument. Consequently, any proof of homogeneity must ensure it does not rely on results obtained using 
an assumption of homogeneity. Checking that the distribution of matter is homogeneous, while still using an FLRW 
model, is at best a consistency check and not a proof. A genuine proof must allow both the geometry and the matter 
to be inhomogeneous. 

Historically, it has been theorized that it should be possible to use astronomical measurements to directly determine 
the metric of the Cosmos. The idea of co-ordinates based on the observer's past null cone was first suggested by Temple 
[1], who called them "optical co-ordinates." McCrea surveyed observational relations in homogeneous models [2] and 
was the first to determine them, as a series approximation, for inhomogeneous models [3]. Kristian and Sachs [4] 
were the first to consider the determination of a general metric for the Universe from observational data. They used 
series expansions near the observer in powers of diameter distance (and density field) for redshift, image distortion, 
number density, and proper motion. They considered the case of dust and estimated the parameters in their series, 
concluding that homogeneity was not proven. Ellis et al. [5] provided a major review of this problem. They showed 
it is not possible to obtain the spacetime metric without assuming Einstein's field equations, and they also raised 
the problem of source evolution. Further papers by Stoeger, Ncl, Maartens, Araujo, and others [6-14] considered 
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the fluid-ray tetrad, the spherically symmetric case and its perturbations, various solutions, and the origin conditions 
in these co-ordinates. Mustapha et al. [15] showed the diameter distance and the density observed on the past null 
cone can be significantly distorted by inhomogeneity. Ribeiro and Stoeger [16] considered the inclusion of a galaxy 
luminosity function, and checked catalogue data for consistency with the theory. A follow-up paper of Albani et al. 
[17] showed that analyses of galaxy statistics are strongly affected by the distance definition used. 

Mustapha, Hellaby, and Ellis [18] considered data that are isotropic about us, and described an algorithm that 
could determine the specific Lemaitre-Tolman-Bondi (LTB) model [19-21] that reproduces the observations. Thus, 
they claimed that any reasonable sets of observations of diameter or luminosity distance against redshift, and number 
counts against redshift, can be fit by an LTB model. They also showed that the effects of source evolution, cosmic 
evolution, and radial inhomogeneity are mixed together on the past null cone, and it is difficult to distinguish them. 
However, Hellaby [22] suggested a method of checking theories of source evolution in inhomogeneous geometries using 
multi-color observations. Celerier [23] considered a series expansion on the past null cone of an LTB model, and showed 
that the supernova data can be reproduced. Recently Celerier reviewed [24, 25] the need for inhomogeneous models 
in relation to observations, concluding that the standard methods of fitting a homogeneous universe to the data are 
open to significant systematic errors, and noting that fitting observations with inhomogeneous models is a significant 
challenge. Bishop and Haines [26] approached the Kristian and Sachs program as a time-reversed characteristic initial 
value problem, but they didn't solve the problem of the maximum in the diameter distance, and they only tested 
their numerical code with Einstein-de Sitter data. 

Nevertheless, other than very limited series expansions, the various proposed methods have never been implemented 
with real data. Therefore we have begun a project to develop a numerical procedure that can do the necessary 
calculations and handle observational data. Thus, when sufficient data become available in the not-too-distant future, 
we will have a numerical routine capable of determining the Universe's spacetime structure. 

The algorithm suggested by Mustapha, Hellaby, and Ellis [18] was implemented as a numerical procedure by Lu 
and Hellaby [27, 28]. The process of turning the algorithm into code raised some important issues that had not 
been obvious when discussing the theory. These include the fact that the real data are discrete, there are no data 
at the origin itself, and the differential equations to be solved are singular at the maximum in the areal radius, as 
well as the problem of the choice of appropriate numerical methods, given that the differential equations contain 
both observational data and the functions being solved for. The program was tested using ideal data, for a variety 
of homogeneous and inhomogeneous models, and successfully recovered the original models, thus demonstrating the 
viability of the method. 

The main reason for working with spherical symmetry in these initial stages is merely to keep the theory and the 
numerics relatively simple. In the long run, this assumption will be dropped. However, there are three observational 
reasons why spherical symmetry is a good place to start. First, we note that we are at the center of our past null cone, 
so it makes sense to consider spacetime in terms of spherical co-ordinates. Second, the Universe does appear nearly 
isotropic on large scales, but radial homogeneity is not easy to verify because of the finite travel time of light and the 
miniscule duration over which cosmological observations have been made, so it is more urgent to determine the radial 
variation of the metric. Finally, there is no deep all-sky redshift survey at present, and the zone of avoidance is likely 
to be a gap in any survey for the foreseeable future. 

If observations are made over the whole sky and the Universe is truly homogeneous on large scales, then using 
a spherically-symmetric model and binning the data over all angles should smear out any structures and lead to a 
radially-homogeneous metric at radii larger than the largest structures. At smaller radii the spherical shells may 
be enclosed entirely within a void or dominated by a supercluster and would appear inhomogeneous. Also, as a 
first check on isotropy, measurements made in different directions on the sky could be reduced separately, using our 
spherically-symmetric procedure, to look for angular variation. 

The aim of the present paper is to extend the previous numerical routine of Lu and Hellaby [28] to be able to 
handle observational uncertainties. We study the stability of the differential equations used in the numerical routine, 
and how the uncertainties in the data should be propagated through the integration of the differential equations. We 
approximate the effect of random errors in observational data by adding Gaussian deviates to test data generated from 
model universes. Using this simulated data, we develop various impovements to the original numerical procedure, 
especially for calculating the derivatives of the areal radius, for dealing with the regions near the origin, and for 
matching the series expansion used about the maximum in the areal radius with the numerical integration regions on 
either side. Finally, we study the effect of systematic errors on the numerical algorithm and show that it is possible 
to correct for an overall systematic error in the distance scale or in the density (or if both are present, one may be 
corrected relative to the other to obtain a self-consistent result). Hellaby [29] showed that there is a relationship 
between the enclosed gravitational mass M, the areal radius R, and the cosmological constant A at the maximum in 
R that is independent of any inhomogeneity: this relationship is the key to detecting and correcting for systematic 
errors. 

The main concern in this paper is the development of a basic numerical procedure that can handle input data 
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containing statistical and systematic errors. Of course, considerable processing is needed to extract the data functions 
needed by our numerical algorithm from the raw astronomical observations — luminosity functions, K-corrections, 
evolution of galaxies in luminosity and mass, galaxy mergers, etc. — and many sources of error must be considered in 
assessing the uncertainties. Though these must be carefully addressed in the future, for the present we assume that 
we are given data in the appropriate form with known uncertainties. It seems that type la supernovae are establishing 
themselves as very reliable standard candles out to redshifts well above 1, and when a sufficient number are measured, 
they may well provide the luminosity distance data for this project, or the basic calibration for a range of sources. 
Determination of the density of matter in redshift space, which is composed of the number density of sources and the 
mass per source, is much more problematic. The determination of masses from dynamics does extract the desired 
gravitational masses, but requires detailed observations, is not easy to do in large numbers, and can only be performed 
where there is luminous matter. It would be very convenient if there were a well-established tracer of mass that could 
be used instead of trying to count everything. There are some suggestions, such as luminous red galaxies, but it is 
very hard to demonstrate a certain type of source is a reliable tracer until a full mass census is taken over a region 
including several structures of all sizes. If the tracer sources are too widely dispersed, this would limit the scale on 
which inhomogeneity can be detected. Although current cosmological data are well short of the required quality and 
quantity, we expect dramatic improvements in accuracy and completeness in the not-too-distant future. 



II. THE MODEL AND THE NULL CONE 



We use the spherically-symmetric inhomogeneous-dust model (the "LTB" model) discovered by Lemaitre, rediscov- 
ered by Tolman, and studied by Bondi [19-21]. We adhere to the notation used by Lu and Hellaby [28]. The metric 
is 



ds^ 



l\2 



1 + 2E 



(1) 



where dVL^ = d9'^ + sin^ OdcjP, R{t, r) is the areal radius, a prime denotes partial differentiation with respect to r, and 
the free function E{r) > —1/2 is a localized geometry term. From the Einstein field equations we get 
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where R = dR/dt, and 
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where M(r) is a second free function that gives the gravitational mass within a comoving shell of radius r. Here E{r) 
also plays the role of the local energy per unit mass of the dust particles. For the present we take A = 0, postponing 
more general considerations for later work. The solutions of (2), in terms of parameter 77, are 
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(5) 
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for hyperbolic, parabolic, and elliptic evolution respectively. (Near the origin, where E 0, the type of evolution 
is determined by the sign of RE/M or E/M'^/'^.) These solutions contain a third free function tBir), which is the 
time of the big bang locally. By specifying the three free functions — M(r), E{r), and tsir) — an LTB model is fully 
determined. Between them they provide a radial co-ordinate freedom and two physical relationships. 



A. The observables and the differential equations 



The background theory is presented in Section 2 of Lu and Hellaby [28] , and here we merely summarize the essentials. 
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The diameter and luminosity ^ distances are 



R 



D 

J ' 



(7) 



where D{z) and L{z) are the true diameter and absolute luminosity of a source, 6 and i are the corresponding angular 
diameter and apparent luminosity, c?io is 10 parsecs, and they are related by the reciprocity theorem [30] ^ 



(8) 



The redshift-space number density of sources n, measured in number per steradian per unit redshift interval, is related 
to the density p by ([18, 28]) 



£>2- da; 
R p = pn— , 
dr 



(9) 



where fi is the mass per source. 

Light rays arriving at the central observer follow ds^ = = d6^ = d0^, so the past null cone of the observation 
event (t = tQ,r — 0) satisfies 



dt 
dr 



R/_ 

W ' 



W = VT+2E 



(10) 



and we write the solution t = i{r) or r = r{t). We denote a quantity evaluated on the observer's past null cone with a 
hat on top or as a subscript, for example R{t{r), r) = R or [i?] a, though this will often be omitted where it is obvious 
from the context. 

We can use the radial co-ordinate freedom to choose 



dt_ 

dr 



-1 , i.e. R' , 



on the observer's past null cone, so that the solution to (10) is 

i=to-r . 



(11) 



(12) 



Note that (11) and (12) and the following differential equations only hold for the single null cone with apex (tg, 0). 
Since the co-ordinate r is not an observable, all r derivatives are converted to z derivatives, defining 



and using 
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Then from Lu and Hellaby 
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the differential equations are 



dR 

dr dz 
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(15) 
(16) 



^ In Kristian and Sachs [4] a "corrected luminosity distance" was defined to be the same as the diameter distance. Some authors have 
called this the "luminosity distance," which has led to a confusion of terminology and sometimes to incorrect definitions. 

^ This was first shown by Etherington [31]. Kristian and Sachs [4] cite a private communication from R,. Penrose for a general proof of 
the reciprocity theorem, but we are not aware of any publication of that work. 
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and 

deb ( 1 ^ + 

Equations (13), (17), (15), and (16) constitute the differential equations (DEs) to be solved for (p{z), r{z), M{z), 
and E{z). Then tB{z) follows from (4)-(6) and (12). Knowing r{z), M{z), E{z), and tB(z) fully determines the LTB 
metric that reproduces the given R{z) and fi{z)n{z) data. 

The LTB origin conditions were thoroughly examined in Mustapha and Hellaby [32] , and Lu and Hellaby [28] give 
their application in this case. 



B. Apparent horizon 



The locus where the past null cone crosses the apparent horizon is also where the areal radius is maximum, i.e. 
the maximum R is R = Rm at z = Zm- Since AR/dz = here, the DEs (17) and (15) with (16) become singular. 
However, it is obvious from (16) that where dR/Az = we must also have 



Rm ~ 2Af,„ , 

since W is arbitrary (see Krasihski and Hellaby [33], Hellaby [34]). Similarly (13) and (17) show that 



A?R 



dr^ 



A^R 



Az^ 



(18) 



(19) 



since we don't expect Az/Ar or d^r/dz^ to be divergent here in a general LTB model with co-ordinate choice (11). 
Indeed, (18) and (19) are exactly what happens at in the FLRW case. So although there are no divergencies at 
Rm, the numerics break down. In Lu and Hellaby [28] this was overcome by doing a series expansion in Az = 2 — z,„, 
and joining the numerical and series results at some z value Zj, < z^ — see sections 2.6, 3.3, and appendix B of their 
paper. As pointed out by Hellaby [29], this phenomenon is not merely a cosmological curiosity. At this locus, and 
no other, there is a simple relation between the diameter distance do — R and the gravitational mass Mm that is 
independent of any inhomogeneity between the observer and sources at this distance: 



2M„ 



Rn 



ARl 



(20) 



or (18) if A = 0. (However, he redshift Zm at which this occurs is model dependent.) Therefore the maximum in 
R provides a new characterization of our Cosmos — the cosmic mass. More importantly for this paper, it provides a 
cross-check on the numerical integration, since the M value obtained from the numerical integration must agree with 
that deduced from the measured Rm using (18) or (20). For this reason, it plays a prominent role in what follows. 



C. Outline of numerical procedure 



Since the long-term intent is to use observational data as input, some important practical considerations arise, as 
discussed in Lu and Hellaby [28] . The above analysis assumes continuous functions and differential equations, whereas 
real cosmological data consist of a large number of discrete measurements of individual sources (e.g. galaxies); and 
numerical integration involves discrete steps of finite size. 

Measurements of the magnitudes (or luminosities) and masses of individual galaxies involve considerable uncertainty, 
so it is advisable to average over a significant sample, and this requires putting the data into redshift bins. (At present, 
bins of width Sz = 0.001 are being used.) While redshift measurements are relatively accurate, the observed redshifts 
include peculiar velocities, so again bin averages reduce the uncertainty ^. The very- low z bins ^ are strongly affected 



^ Actually, the diameter distance is unaffected by peculiar velocity, even if it is calculated from the luminosity distance and the redshift 
using the reciprocity formula (8). Therefore, if measurements of diameter distance or luminosity distance become sufficiently accurate 
out to large redshifts, it would be preferable to bin the data hy do = R. The main difficulty would be to distinguish sources just closer 
and just farther than the maximum in R. 

* If we could detect all sources, the low-z bins would have very little data, but in fact the data are far more complete at low z, so the 
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by peculiar velocities, and there are no data at the origin itself, which makes it difficult to set the initial conditions 
for the DEs. Thus, we have to estimate the origin values from the data in the low-z data bins and a series expansion 
of the model. Therefore the procedure for integrating down the PNC has four main stages: the origin fitting, the first 
numerical integration interval, the near maximum series, and the second numerical integration. 

Because the DEs involve the first and second derivatives of the observational quantity it is evident that noise 

in this function could create very large fluctuations in the derivatives. Thus, in addition to binning the data, it may 
be necessary to smooth out statistical fluctuations by fitting a curve to the binned data. 

Each of the procedures — binning, smoothing, numerical integration procedure, etc. — is a potential source of unin- 
tended bias and obviously determines the smallest scale on which inhomogeneity can be detected. So it is likely that 
our thinking about the best way to handle all these issues will evolve as the project proceeds. 

The next two sections consider the effect of statistical uncertainty in the data functions R{z) and 47r/i7i(z). This 
affects both the uncertainty in the output functions M(z), E{z) and tsiz), and the stability of the solution procedure. 



III. STABILITY OF THE DIFFERENTIAL EQUATIONS 
A. Equations for r and 

The co-ordinate distance r is calculated by integrating (13), once = dr/dz has been determined from (17), 

d^ 11 Rzz 4:Trfin(f) \ 



(21) 



dz yi + z 4 J 

where z subscripts denote differentiation with respect to z. The expression for the error in (j)^ due to the error in (j) is 

A^. = (I + ^) (22) 

which holds only for Acf) <C (j). At small z, the 4:TTfin(j)/{RRz) term is insignificant relative to the other two terms in 
4>z at z = 0, and Rz > before the maximum in R, so this leads to the integrated values of (j) being stable so long 
as (j}z is negative and (f> is positive. When (f>z is stable, if is over (or under) estimated, then the factor of (f> in (j)z 
leads to (pz being under (or over) estimated so that </> opposingly decreases (or increases) toward the correct value. 
It is possible for to become unstable if (j)z becomes positive; however, this should only occur if there are strong 
inhomogcneities, and only over a small range of z, meaning the instability can only grow over a short range before 
becomes stable again and corrects itself. 

Beyond the overall influence of the (j) factor in 0^ , there is an additional in the third term of 0^ that can lead to 
instability before the maximum in R. As the maximum Rm is approached, the 4,1: jjmcj) I RRz term grows faster than 
0^/0, and after the point where 



Rzz ( 1 STTfincj) 



Rz \l + z RR^ 



(23) 



the DE becomes increasingly unstable towards Rm- Near the maximum, the terms in (21) involving 1/Rz both 
diverge, but they ought to cancel with ideal data and calculations, because Rzz = 47r/^n0/i? at Rm- Since the sign of 
Rz must change from positive to negative at i?„i, while all other quantities maintain the same sign, then returns to 
being stable after Rm- 

Since the values of r are integrated using the values of 0, the stability in r depends on the stability in 0. The main 
difference is that if diverges in one direction and later stabilizes back to the proper value, there is a cumulative 
impact on the integrated values of r, so the consequence of any instability in is potentially worse for r. The 
integration relies on the use of a series expansion to get past the maximum in R- Provided the integration switches 
over to the series expansion well enough in advance of the maximum, then should not diverge significantly from the 
correct value, and the cumulative effect on r should not be significant. 



opposite is true. 
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B. Equations for M and E 



The differential equation (15) with (16) for M is given by 

az 

where 



, R. ( 2M\ 6 

H \ R 2R, 



and so the stabihty equation for AM <C M is 




R J 2R,\ RR, 



AM 



(24) 



(25) 



(26) 



Since W contains a term that goes as —M(j)/{RRz), if M is over (or under) estimated, then W, E, M^, and M 
are opposingly under (or over) estimated when R^ is positive, making the determination of M from W stable and 
the determination of W (and E) from M stable. Barring the presence of strong inhomogcncitics, R^ should remain 
positive out to the maximum in i?, so M and W should remain stable with each other up to the maximum. If there 
are strong inhomogeneities, these should only make Rz go negative over short regions of the integration, meaning 
the instability can only grow briefly before M and W become stable and return to their proper values. The values 
of W, E, Mz, and M also depend on the integrated values of 0, so if (j) becomes unstable, as it may do before the 
maximum, then W and M can potentially become unstable due to but as discussed previously, the series expansion 
of (f> around the maximum should bridge this instability in </). 

After the maximum in R, Rz should remain negative (except in the presence of strong inhomogeneities), which 
means over (or under) estimates in M lead to over (or under) estimates in W and E, and over (or under) estimates in 
W lead to over (or under) estimates in Mz and M. Thus, M and W feed back on each other and lead to instability. 
This suggests an alternative method is needed to determine M and E after the maximum in R is reached. Assuming 
the data do not extend significantly beyond the maximum, one method is to create an extended series expansion in 
W based on the Anfin and R data from the maximum to the outer boundary, and then integrate values of M based 
on the W series expansion so that M and W do not feed back on each other. If the data extend significantly beyond 
the maximum, and no alternative method is found, we must expect very large uncertainties in W and M. 



IV. NUMERICAL HANDLING OF STATISTICAL ERRORS 
A. Simulating errors with Gaussian deviates 

In order to make sure the code can handle real data, statistical errors must be added to the distance determinations 
and galaxy number counts for the test data. 

The actual number of sources N per bin is equal to AnnSz (where 6z is the redshift interval of the bin), since n is 
the number of sources per steradian and per unit redshift interval. From the test universes, the number of sources 
is calculated as 

N^—t- , 27 

M 

where /i is the assumed mass of a typical galaxy — 10^^ Mq. The program then calculates the expected errors in the 
R and AiTfin values by calculating the l-cr errors and using the elimination method to add Gaussian deviates (e.g. 
[35]) to the data in a ±5cr range. 

The random error in a single R measurement is assumed to be of order 10% so that the random error AR in the 
R value for a bin of N galaxy sources is 



(28) 



8 



(Although current luminosity-distance measurements are mostly not this good, it is expected the accuracy will improve 
considerably in the coming years.) 

The random error in the A-nfin measurements is 



A(47r/xn) = — JiApiNY + (AN ^if = — V(AmAiV)2 + (AAf/x)2, (29) 
oz dz 

where AA^ = ^/N, and Aji = Am/ A A relates the error in the mean mass per source /i and the error in a single mass 
measurement m. The uncertainty in A^ is calculated using the Poisson error, since if the galaxies occur randomly, 
then the uncertainty in the number that should be observed goes as Va". (In reality galaxies are clustered, but since 
we are binning over thin spherical shells, clusters will be averaged out at large radii.) Assuming Am is less than (or 
of the same order as) then the random error in the A-n^n measurements is approximately 

A(4.;.n) = ^ = (30) 



B. Calculating derivatives of R 

With errors in the R values, it is no longer possible to use the original method of Lu and Hellaby [28] of calculating 
the first and second derivatives, which were obtained from the first and second differences in R between adjacent 
redshift bins. Due to errors in the data, the noise in the second derivative ends up being of order 10000%. This causes 
the integrated values of </> to vary rapidly up and down and quickly diverge shortly after the instability in (p develops 
(before the maximum in R). 

Thus, in order to reduce the noise in the derivatives, it is necessary to calculate least-squares fits to the R data. 
Unfortunately, since the functional form of R{z) cannot be assumed, simply fitting a single polynomial is unlikely to 
yield a good fit over the whole range of data. Rather, the values of the first and second derivatives that best reflect 
the local data are obtained by fitting a polynomial to the 2fc -I- 1 redshift bins centered on each z,;, covering the range 
Zi ± k5z. Since variations in the second derivative should exist, especially if strong inhomogencities are present, it is 
necessary to use at least a cubic polynomial, so that the second derivative is allowed to vary over the range of the fit. 
Near the origin and outer boundary, it is not possible to have a range of data that is centered exactly on the point in 
interest, making it even more important to have at least a quartic fit so that the unequal upper and lower ranges of 
data do not bias the calculation of the second derivative. 

When there are strong variations in Rzz, it is necessary to either reduce the range of data used in the fit or increase 
the order of the fit so that they do not get smoothed out. In either case this leads to more noise in the fit. The ideal 
data range and polynomial order varies depending on the test universe. If there are no strong variations it is better 
to use a lower order polynomial (such as a quartic) and a large range of data (up to 300 redshift bins centered on 
the current z value) to reduce the noise, but when there are strong variations, it is better to increase the order or 
reduce the range to allow some noise, in the interest of fitting the actual variation. It is obvious if real variations 
are present, since these are noticeable even in the lower-order fits where the noise is negligable. Exactly how much 
the range should be reduced or the order of the polynomial should be increased can be ascertained by changing the 
fit until the real variations cease to become stronger or when the noise starts to become a significant fraction of the 
magnitude of the real variation. In practice, for the test model with strongest inhomogeneity, the range could be as 
small as 30 5z bins and the fit as high as 7th order. 

More precisely, the best polynomial fit to the data can be ascertained by comparing the ratio of with the number 
of degrees of freedom (NDF) , where is the sum of the squares of the R residuals divided by the 1-a error of each R 
value, and the NDF is the number of data points minus the number of polynomial coefficients being fit. The expected 
value of is the NDF, so the ratio of the two that is closest to 1 yields the best fit to the data. A ratio greater than 
1 is not adequately fitting the data, while a ratio smaller than 1 is fitting the data too well, which means it is likely 
fitting the noise and is not the most realistic fit to the data. 

It should be noted that the least-squares fitting of R versus z minimizes the error in the R values, so it assumes no 
error in z. Any actual errors in z would bias the fit. While the observational errors in z should be small, the peculiar 
velocity components of the observed redshifts may lead to errors. 

Another consequence of this smoothing is that random fluctuations in the R data create quasi-periodic variations 
in the smoothed versions of i?, i?^, and Rzzi which have a period determined by the z range chosen for the fit. These 
manifest themselves in (j) as oscillations with an amplitude that diverges as the maximum Rm is approached. 
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C. Using an origin LTB series expansion 



The integration of the differential equations must proceed from an initial value, which in the original code of Lu 
and Hellaby [28] involved using an FLRW series expansion at the origin to determine the derivatives in the first two 
redshift bins. Since a regular LTB model is FLRW-like at the origin, this is not a bad first approximation. However, 
the fractional errors in R and AiTfin are largest near the origin, due to the smaller number counts that occur in a bin of 
given thickness Sz, and also, the R derivatives are more susceptible to error near the origin. The determination of the 
bang time is more difficult near the origin, because this depends on 2ER/M , while E, R, and M are all approaching 
zero and highly susceptible to numerical errors. Thus, it is preferable to use the series expansion for more than just 
the ffi'st two bins, which necessitates the use of an LTB series so that the region around the origin is not constrained 
to be FLRW. Approximating R and Anfin near the origin by 5th order power series in z, the series coefficients are 
found by a least-squares fit to the data in the first 50 bins (which extends to z = 0.05 — a distance scale at which there 
should be adequate data to begin integrating from bin to bin without large numerical errors). The functions cj), M, 
and E are also expressed as power series, and their coefficients are derived from the DEs. The values of the functions 
at the end of the origin series then provide the initial conditions for the numerical integration. 

Series expansions of the R and At: fin data near the origin in the form 

R{z) ^ Riz + R2Z^ + R^z^ + RiZ^ + R^z^ (31) 



and 

4:TTfin{z) = K2Z^ + Ksz^ + K^z^ + K^z^ (32) 
are used to calculate the series expansions of r, </>, M, and 2E from (17), (13), (15), and (16), as 

r{z)=R,z+{^^ + R2)z'- + {^ + R, + !^)z^+[^+^-^ + R, + % + i^)z' 
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and 
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+ V^8i?r + tr^ ^ iRP" J + {2qr: + 5i?r^ ) sr: 

In Figure 1 , the relative error in the bang time near the origin is shown for both the previous method and the new 
method. The large numerical error that occurs with the previous method is greatly reduced by using the LTB series 
expansion. 



D. Bridging the maximum in the areal distance 

Since, as noted in Section II B, R^ goes to zero at Rm (the maximum in R), and the DEs for (f> (17) and M (15) 
and (16) both depend on 1/Rz, it is necessary to use a series expansion in the neighborhood of the maximum. Also, 
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FIG. 1: Comparison of the relative errors in the bang time ts versus z using the old (crosses) and new (asterisks) methods. 
The new method fits an LTB series expansion to the first 50 points, before starting the numerical integration, thus allowing 
the solution to be inhomogeneous while avoiding the numerical difficulties near the origin. This plot is for an inhomogeneous 
universe with origin Hubble and deceleration parameters ho = 0.72 and go = 0.6, and with JIa = 0. 



the (j) DE becomes unstable before the maximum, so this series expansion is necessary over a range of more than just 
a few points. 

Appendix A gives the series solutions for 0, Af , and W in powers of Az ~ Zm. — z. The coefficients in the series are 
fully determined by the DEs, but the M and W coefficients are interrelated and contain one undetermined coefficient, 
the linear term in the M series, or equivalently, the constant term in the W series. In fact every coefficient of the M 
and W series depends linearly on this factor, so its value is easily fixed by matching these two series to their numerical 
values at some cross-over point z = Za < z„i, where the series and numerical (/) curves are in closest agreement, as was 
done by Lu and Hellaby [28]. A second matching is done at z = zj > Zm, where the near- maximum series ends and 
the numerical integration is re-started. 

When statistical errors are present however, the oscillations in (f> that grow as is approached make the matching 
of the numerical and series values problematic, since they will tend to intersect at multiple points. Thus, the code 
has been modified such that the two <f> curves are matched by calculating the absolute value of the differences between 
the curves over a range of the data (the same range used in the least-squares fit of R) and matching them at the z 
value where the total difference is minimized. This generally leads to a matching near where the instability sets in, 
thus excising the region of diverging (jj error. 

In the original code of Lu and Hellaby [28], M was matched at the start of the series, and W was matched at the 
end of the series. While M and W are stable with each other before Rm, the fiuctuations in influence W, which 
opposingly influences M, so typically there is a small error in M at the first matching point. This M matching error 
biases the slope of the M series, and in turn the constant term of the W series. Thus, there is generally a small jump 
in W at the first matching point and a small jump in M at the second matching point, which in turn leads to a rapid 
adjustment in W once the numerical integration re-starts. 

However, it turns out that it is possible to find values of M and W that are consistent with each other at the first 
matching point so that the small error in M does not lead to overcorrections. The mean of the integrated value of 
M and the value of M that would result from matching the integrated value of W is used to determine the series- 
connection value of M. Similarly, the mean of the integrated value of W and the value of W that would result from 
matching the integrated value of M determines the series-connection value of W. These new series-connection values 
of M and W are consistent with each other in that calculating the M-series matching using the mean W value for the 
matching leads exactly to the mean M value, and vice versa for the mean W value. Also, these values are stable in 
that calculating the matching with "numerical" values above (or below) the mean results in series-connection values 
closer to the mean values. This can be explained by the fact that M and W are stable with each other before Rm, so 
there should be stable series expansion values of M and W that are consistent with each other before Rm- 

Plots of the relative error in M and W in the region containing the matching before Rm appear in Figure 2. Using 
this new matching method reduces the jump in W at the first matching so that it does not overcorrect, and it leads 
to more accurate values of W throughout the series expansion. It also leads to more accurate values of M in the 
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FIG. 2: Comparison of the relative errors in W (left) and M (right) versus z in the region before Rm, using the old (small 
crosses) and new (large diamonds) methods for matching from the numerical values to the near-maximum series expansion 
values. These plots are for an FLRW universe with ho — 0.72, go = 0.22, and = 0, with statistical errors added to the data 
to simulate observational errors. Simply matching M can lead to a large jump in W and a divergence in M, whereas the new 
method leads to small jumps in both M and W that always bring them toward the correct values. 



near-maximum series, and much better values of M at the second matching. This is an important improvement over 
the original code, especially when there are errors in the data that can lead to larger jumps in M and W. 



E. Propagation of uncertainties 



To estimate the uncertainties, the observational errors need to be taken into account in how they propagate in the 
integrations. Where a series fit to the data is performed, the estimated errors come from the covariance matrix of the 
fit. The uncertainties in the numerical integration of r^, 0^, W, and Mz are 

Ar, = A0, (37) 
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(40) 



assuming that the errors arc uncorrclated with each other so that they are added in quadrature. 

The (j)z values are used to calculate each value from the previous (j) value, but how the error in (f> accumulates 
depends on the stability in (j). It would be expected that when (j) is stable this error would not grow from step to step, 
while when (j) is unstable the errors from all the unstable steps would need to be added together directly. In practice, 
calculating the uncertainty is more complicated, due to the smoothing of the R data. The errors in the smoothed 
Rzz values tend to fluctuate up and down with a period = 2k6z equal to the range of data that is used to perform 
the least-squares fit. These fluctuations in the R^z errors induce fluctuations in the (j) values, so even when (j) is stable 
it tends to oscillate with the same period as the range of the least-squares flt. However, over a range > the errors 
due to the fluctuations tend to cancel, and so should not be added continuously. 
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FIG. 3: Relative <j) error versus z, with the estimated error (large points) plotted alongside the actual error (small points) for 
comparison. Gaussian deviates are added to the data, and the estimated uncertainty calculated, along with the actual error 
from the ideal data. It can be seen that the error in (f) grows before the maximum as it becomes increasingly unstable, and then 
decreases after the maximum as it becomes stable again. This plot is for an inhomogeneous universe with origin parameters 
ho = 0.72, qo = 0.6, and Qa = 0. 

The actual calculation of the cj) uncertainty proceeds by taking the errors of all the previous steps that may 
potentially be continuously fluctuating in one direction (starting either from the initial point of integration or at a 
point that is previous to the point of interest Zi by half the range of the R fit, C/2) and adding these in quadrature 
(since the errors at these steps are unlikely to all be of the same sign) to get the error due to the fluctuations At/)/: 



and combining this error with the error A0s from the end of the series leading up to the numerical integration by 
adding them in quadrature: 



The error A0 / is due to all the uncertainties in the data (38) over half the smoothing range (or less if the start of 
the numerics occurs within that range). The junction of the numerical region to the near-maximum series expansion 
generally occurs just where (f> becomes unstable and A(f> starts to diverge, so it is essentially unnecessary to ever start 
adding the full A0 errors directly instead of in quadrature. 

It is assumed that the errors coming out of the origin series are zero, since in practice the actual errors are always 
miniscule compared with the estimated error beyond the first few points of the series. Beyond the maximum, in the 
second numerical integration region, the error due to the near-maximum series is not negligible, but this is slowly 
corrected by the natural stability of the (j) DE, so that it decays from step to step as 



where A(/)[io] is the error at the end of the near-maximum series, and A0* is calculated from (22). reflecting how the 
error in acts to stabilize itself. Figure 3 contains a plot showing how the estimated error compares with the actual 
error in the propagation of the uncertainties. 

The integrated r values come directly from the values, but the uncertainties have to take into account the 
uncertainties in all the (f> values. Even when cj) is stable and self-corrects when an error emerges, the integrated values 
of r all continue to be off by the amount gained or lost while the values of </> were in error. When cj) fluctuates above 
and below the correct value, then the errors in the r values grow like a random walk. Thus, the error in r is calculated 
by adding all the previous errors due to in quadrature to account for how far it should have wandered due to the 
random high and low values of (j). 




(41) 




(42) 




(43) 
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The integrated values of M are stable before the maximum in R and then become unstable after the maximum due 
to the feedback with W . Thus, before the maximum, it is assumed that there is no M error in the calculation of the 
W error, and no W error in the calculation of the M error, so that these errors do not feed back on each other and 
depend solely on the An^n, 0, i?, and errors. After the maximum, the W errors are included in the calculation of 
the Mz errors, and if the W values are being determined from the integrated values of M (rather than an extended 
series expansion for W), then the M errors need to be included in the calculation of the W errors as well. 

V. DEALING WITH SYSTEMATIC ERRORS 

It is very likely that there will be systematic errors in the observations as well as the evolution functions. We here 
investigate the effect of such errors, and show that the data functions R{z) and 4:TTiJ,{z)n{z) must satisfy a consistency 
condition. We demonstrate how a discrepancy in this condition can be corrected in order to make the data and 
the numerical result self-consistent. However, it is not possible to determine from this process how much error is 
attributable to R and how much to Airfin^ or to discern whether the systematic errors vary with z. 

A. Systematic error in the distance scale 

We first consider a constant percentage error in R only. In this case, the origin values of (f> start correspondingly 
high or low, since 4> goes as R/z near the origin. Then by (21) the numerical integration of (j) continues to be off since 
(j)z depends on a Rzz4'/Rz term. When the At: ^ncj) / [RzR) term in (f>z becomes significant, (p starts to correct itself 
since the errors in l/(RzR) bias the integration of in the opposite sense, and the integration may even overcorrect 
if it proceeds long enough. 

By (A3)-(A7) the erroneous R values produce correspondingly high or low values in the near-maximum series for 
(j>, but because of the correction that takes place before the maximum, there is generally a discrepancy between the 
numerical and series (f) values so a close matching is not possible. This discrepancy is a clear indicator of an overall 
systematic error. Beyond the near-maximum series, the numerical integration is stable and starts to correct the value 
of (f>, leading to a discontinuity in the slope of the <j> curve, providing a secondary indication of a systematic error. 
Since the integrated values of r depend on the (f> values, then if (f) is too high (or low) over most of the run, then r 
tends to be correspondingly high (or low) over the whole run by approximately the percentage of the error in the R 
values, but it is not obvious from the r{z) curve that a systematic error is present. 

The presence of a jump in (f) at the start of the near-maximum series not only flags a systematic error, but also 
provides the means of correcting it. Thus, simply by systematically increasing or decreasing all of the R values until 
the jump in (j) is minimized, the best correction for R can be found. The adjustment of the R data necessitates a 
repeat of the entire integration from the origin to the maximum series matching point. The effects of a systematic 
error in i? on (/) and r can be seen in Figure 4. 

A systematic error in R is also apparent in W. Initially, the values of W and M are unaffected by the R values. The 
values of (j), R, and Rz are generally all correspondingly high or low, so when the —M(p/{RRz) term in (25) starts to 
become significant, W becomes correspondingly too high or too low (and M in turn also becomes slightly too high or 
too low). However, when matching the numerical integration to the near-maximum series, the change in the value of 
i?m/2 is greater than the percentage change in M , which influences the slope of the M series and the constant term 
of the W series such that M and W end up being correspondingly high or low in the ensuing numerical integration. 
Since the calculation of depends on E and M, the jumps in W and M also typically lead to jumps in tB- The 
effects of a systematic error in R on W, M, and ts can also be seen in Figure 4. 

It should be noted however that only a constant correction can be made this way, and a z-dependent systematic 
error cannot be distinguished from a constant one. The ability to detect an accumulated systematic error and correct 
for it is due to the unique geometric properties of Rm noted by HcUaby [29] , and this locus is an exception to the 
theorem of Mustapha et al. [18]. 

If the systematic percentage error in R grows with z, the effect at low z is qualitatively similar, but less pronounced 
near the origin. However, the maximum in R is changed not only in magnitude Rm, but also in redshift Zm- If this 
occurs, then the series is calculated opposingly low (or high) due to the higher (or lower) value of Anfin at Zm- 
Then the subsequent numerical integration only asymptotically corrects the value of </>. Thus, the overall influence 
with cj) being successively too high and too low for parts of the integration is that r may end up approximately correct 
at the end of the integration. Both M and W are not much changed before the series matching, but the shift in 
Zm opposingly decreases (or increases) the linear slope of the M series curve more than the high (or low) value in 
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FIG. 4: Systematic errors in the data lead to an inconsistency when attempting to join numerical and near-maximum series 
regions, leading to jumps in 0, M and is- This is illustrated for a systematic error in R only. The functions <f) (top left), r 
(top right), W (middle left), M (middle right), and ts (bottom) versus z are plotted with the exact data (thick line) and with 
various corrections of the systematic error (thin lines). Each subsequent curve represents a 5% difference in the correction of 
the R data. The plots are for an FLRW universe with ho = 0.72, qo = 0.22, and Q.a = 0. 
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Mm = Rm/'2' decreases it. This makes the values of M and W too low (or high) when the numerical integration 
resumes after the near-maximum series expansion. 

B. Systematic error in the redshift space density 

With a fixed percentage error in /in, 4> is initially unaffected, but eventually becomes correspondingly high (or low) 
when the A-k jjincj) / {RRz) term in (21) becomes significant. The near-maximum series for </> is opposingly low (or high) 
at Rm by (A3)-(A7). Coming out of the series, the subsequent numerical values of (j) only asymptotically approach 
the correct value. Thus, with (f) being both too high and too low in turns, the integrated values of r are not as far off 
by the end of the integration as with a systematic error in the distance scale. These effects can be seen in Figure 5. 

The value of M is initially correspondingly high (or low), since it depends on the At: fin values, but W correspondingly 
increases (or decreases) because of the —M(f)l{RRz) term. This causes M to be only marginally high (or low) at the 
matching, but enough that the linear slope of the M series (and constant term of the W series) is opposingly low (or 
high) so that this then continues throughout the final numerical integration. The effects on the W and M values also 
affects the ts values. These effects can also be seen in Figure 5. 

Again, the presence of a vertical jump in (p at the matching from numerical to near-maximum series regions is a 
distinctive feature of a systematic error in the data, which may be corrected by adjusting the fin data until the jump 
vanishes. 

The results are largely the same for an error in fin that systematically grows with redshift, except that obviously 
the effects are smaller if there is no error to begin with and the effects only begin to manifest themselves at higher 
redshift. 



C. Systematic errors in both distance and density 

Realistically, we expect both R and fin to suffer from systematic errors. The plots in Figure 6 show the effects 
when the R values are all 10% too large and the Anfin values are all 10% too small. If a correction is made to the R 
values only, this leads to an overcorrection such that r and M are 10% too small and is 10% too large, compared 
with the original model the data came from. If a correction is made to the Airfin values only, this makes r and M 
10% too large and 10% too small. Once one correction has been made, there are no jumps left in the data, so it 
is not apparent any further correction needs to be made. 

With identical constant percentage errors (say p) in both R and Anfin, (p starts out correspondingly high or low, 
the (j) series is only slightly off, and then the (j) integration asymptotically approaches the correct value at the end. 
Thus, r is too low by nearly the same percentage error, over the whole run. The low Anfin values tend to make 
M low, and the errors in the —M(j)/{RRz) term for W essentially cancel out so that W is correct. Since both the 
M matching value and the M series value at the maximum are correspondingly high or low, then the slope in the M 
series (and hence the intercept in the W series) is approximately correct, so M ends up being off by about p percent 
over the whole run. Thus, only in the special case where both the errors are systematically off by the same percentage 
in the same direction is the simple result obtained that r and M are both off by the same percentage. In that case, 
since there is no apparent jump in (j), there is no way of knowing that an error correction is necessary. 

Therefore, in making corrections for systematic errors, it is only possible to ensure that the R and fin data are 
consistent at the maximum, i.e. that the accumulated mass derived from the integration agrees with that deduced 
from Rm- 

VI. CONCLUSION 

The idea of extracting metric information from cosmological observations has been investigated theoretically in 
many papers, but until now had never been turned into a practical procedure. Such a project began with Lu and 
Hellaby [28] in which the theoretical algorithm of Mustapha et al. [18] was turned into a numerical program that uses 
galaxy redshifts, luminosity or diameter distances, and density in redshift space, to determine the geometry of the 
Cosmos; that is, it extracts information about the spacetime metric. The program was tested with fake data for which 
the true spacetime metric was known, and successfully demonstrated the basic viability of the procedure. At this 
early stage of development, spherical symmetry is assumed in order to focus on issues such as the development of an 
appropriate numerical method. In developing this procedure, it became clear that there are four distinct integration 
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FIG. 5: Similar to Figure 4, but with a systematic error added to the Aivfin data instead of the R data. The jumps due to 
various corrections of the systematic error in 47r//n are illustrated for the functions (j) (top left), r (top right), W (middle left), 
M (middle right), and ts (bottom) versus z, with the exact data (thick line) for comparison. Each subsequent curve (thin 
lines) represents a 5% difference in the correction of the Anfxn data. The plots are for an FLRW universe with ho = 0.72, 
go = 0.22, and = 0. 
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FIG. 6; Corrected curves with systematic errors added to both the R and 47r/in data (+10% and -10% respectively). The 
plots show the exact data (dashes) and the corrected curves obtained by matching the curves for the R data (dots) or Att^u 
data (solid curve) to eliminate jumps. It is apparent the jumps may be removed by correcting either the R or the An fin data. 
It is not possible to determine which correction should be made, but either result is much more self-consistent than with no 
correction. This is illustrated for the functions (top left), r (top right), W (middle left), M (middle right), and ts (bottom) 
versus z. The plots of r, M, and Ib clearly end up overcorrected by 10% when a correction is made to one variable. These 
plots are for an inhomogeneous universe with origin parameters /lo = 0.72, go = 0.6, and JIa = 0: a strong inhomogeneity is 
present around z — 0.2 : 0.3. 
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regions that require different treatment: setting up the initial values at the origin, a first numerical integration region, 
the neighbourhood of the maximum in the diameter distance, and a second numerical region. 

The original code of Lu and Hellaby [28] has been modified to deal with statistical and systematic errors, and 
we have studied three consequences: how to estimate uncertainties in the metric functions, how data errors affect 
the stability of the integration process, and how a key consistency condition can be used to identify and partly 
correct for systematic errors. Simulated data were generated by adding Gaussian deviates to ideal data based on 
various homogeneous and inhomogeneous model universes, and this was used to test the new program. Significant 
improvements were needed in several aspects of the procedure. Because the DEs to be solved involve Rzz, the second 
derivative of the data function R, it is essential to smooth out observational noise by fitting a polynomial to R over 
a range of z about each point. Without this, calculated Rzz values would suffer from gigantic errors. The estimation 
of the initial values at the origin was improved by using a small-z series expansion of the LTB metric, and fitting to a 
range of redshift bins near the origin. Also, the junctions between the numerical and near-maximum series expansion 
values for M and W = y/TT2E were improved so that numerical discrepancies were minimized. 

An important result is the theoretical and numerical investigation of stability. The problem of integrating through 
the maximum in R, where the DEs are singular, was already solved by Lu and Hellaby [28] by using a series expansion 
about Rrn- The numerical integration of (t>{z) was here shown to be stable away from Rm, and the region of serious 
numerical instability was avoided by the use of this series. The numerical integration of M and W, however, is only 
stable before Rm is reached, and beyond this point errors in M and W feed back on each other. Because of this, it 
will be difficult to extend this integration far beyond unless another method can be found. We applied the quick 
fix of extending the near-maximum series for W to larger z values and combining it with the numerical integration of 
A'/, which removed the instability and gave acceptable results for the test models considered up to z = 3. However, 
such an approach must necessarily become increasingly inaccurate with z. We suspect however that this instability 
is not confined to the LTB model or spherical symmetry, but is a significant issue for analyzing high-z data, though 
the assumption of a homogeneous Robertson- Walker metric may hide the problem. 

Particularly interesting was the effect of systematic errors. There is a relationship (18) (or its generalization (20)) 
between the mass and the diameter distance, that must hold at the maximum in R. With data extending out past 
this maximum, it has been shown that it is possible to correct for an overall systematic error in either the distance 
ladder or the redshift matter density, making them mutually consistent at the maximum. Since a correction may be 
applied to either of the two data functions, R{z) and ^(z)n(z), there is a range of possible results. But the most likely 
systematic error will be in the redshift matter density, because the mass per source will be uncertain due to the fact 
that this mass has to include a component for all the material in the Universe that never collapsed to form objects, 
and for each source that we do see, there will probably be many smaller sources that we do not. Thus, it seems most 
logical to assume the distance ladder is correct and make the correction in the redshift matter density. 

Our broad conclusion is that when sufficient cosmological data become available out towards a redshift z ~ 2, it 
should be possible to calculate the metric for the Cosmos, and thereby check homogeneity, at least within this range. 
With upcoming surveys, such data should be available in the not-too-distant future. At larger redshift values, well 
beyond Rm at roughly z ^ 1.5, there is a difhculty with stability, potentially resulting in large errors, which we believe 
may be a general problem for analysis of high-z data, making it difficult to properly verify homogeneity on the largest 
scales. It is possible that another method might be devised to avoid this instability, but we have not found one. 
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APPENDIX A: THE NEAR-MAXIMUM SERIES FOR THE PNC IN LTB MODELS 



This appendix provides the solution for the LTB arbitrary functions as series in Az = z — z„i near z = Zm where 
the maximum i?,„ = R{zm) occurs. We write out all functions as power series in Az: 



R = Rm + Y. ' 4^^n ^Km + Yl ' f = rm + Y^ r,Az^ , (Al) 

1=2 i=l i=l 

oo oo 

M = Mm + ^ AfjAz* , ^/\ + 2E = W = Wm+Y^ W^I\z' . (A2) 



4=1 
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The coefficients in the R and AiTfin series are found from polynomial fits to the data near the maximum in R, and 
the coefficients of the series for f, M, and W are obtained by substituting these series into the DEs (13), (17), (15), 
and (16). The following extends the results of Lu and Hellaby [28]. We find 

00 = r-i = , (A3) 

rn 

= -2 = f I ^ - ^ - 3i?3'l ^ , (A4) 

02 = fS 

(A5) 

( { A'., 2K^K2 




12X^(1 + Z„,) 4A%„(1 + Zra? 4(1 + Z„ 

Ki 3 



f K2 


^ 4 




4a:2 




1 




1 + Zm 




1 


\ 6A™ 


2(1 + 


2A'4 




5i^™ 





K„,{l + z^) 4(l + z™)2 

i(2i?4) -5i?5 



R2 3R2R3 \ Rn 



>R3 



2Kl K2KI Kl 2K3 



9KI 2Kf^ ml 5A„(l + z™) 
llA'iATa Kl Kl 



2 



18if^(l + Z^) 4A'3,(1 + Z„0 9/^3,(1 + Zra) 



6A'„,(1+Z„)2 l2Kra{l+Z^f 8(l + z„)4 

[ 3JC3 KiK2 3A'f ZK2 
\ AK^ Kl + 8A-3^ + 4if™(l + z„) 

5Kf 



8ir^(l + z„0 
f 4A'2 




20 A „ 20(1 + z„) 



3i?| 26R2R4 

4:Rm 15i?r,i 




(A6) 



(A7) 
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